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Abstract Neutron matter is an intriguing nuclear system with multiple connections to other 
areas of physics. Considerable progress has been made over the last two decades in exploring 
the properties of pure neutron fluids. Here we begin by reviewing work done to explore the 
behavior of very low density neutron matter, which forms a strongly paired superfluid and is 
thus similar to cold Fermi atoms, though at energy scales differing by many orders of magnitude. 
We then increase the density, discussing work that ties the study of neutron matter with the 
determination of the properties of neutron-rich nuclei and neutron-star crusts. After this, we 
review the impact neutron matter at even higher densities has on the mass-radius relation of 
neutron stars, thereby making contact with astrophysical observations. 
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1 INTRODUCTION 


The properties of neutron matter have long been recognized as critical to the 
properties of neutron-rich nuclei and neutron stars. Low-density neutron matter 
critically impacts our understanding of neutron-rich nuclei. Similarly the equa¬ 
tion of state of high-density low-temperature matter is critical in determining 
the properties of neutron stars, including the nrass-radius relationship and the 
neutron-star maximum mass. 

During the past few years interest in neutron matter have been resurgent for 
several reasons. At very low densities, neutron matter is very similar to cold 
Fermi atoms near unitarity (infinite scattering length), since two free neutrons 
are very nearly bound. This enables stringent tests of theories of fermions in 
this strongly-interacting regime ([!]). The equation of state (2,3) and pairing gap 
of unitary fermions (JiJj^, 6) have been calculated precisely and measured very 
accurately. These experiments provided severe tests of the theories; some of the 
calculations (including those reviewed here) provided excellent predictions while 
others were less successful. 

There has also been a resurgence of interest in the properties of neutron-rich 
nuclei. New facilities, such as FRIB, will probe the properties of many hitherto 
inaccessible nuclei. Among the several expected discoveries, the history of nucle¬ 
osynthesis could be mentioned: we still do not know where the heavy elements 
were created. 0 Such experiments will probe the extreme isospin-imbalance 
limit, thereby informing nuclear models and ab initio approaches. Typically den¬ 
sity functionals (8) are used to predict the properties of these nuclei; while some 
components of the density functional are very tightly constrained from stable 
nuclei the extreme isospin limit is less constrained. Many-body calculations of 
inhomogeneous neutron matter can provide further constraints on the parameters 
entering density functionals. This is also important for the inner crust of neutron 
stars, where inhomogeneous neutron matter fills the space between the lattice of 
heavy neutron-rich nuclei. 

Finally, neutron stars have become even more important in recent years. The 
discovery of the first two-solar-mass neutron stars (9,10) provided critical con¬ 
straints on the dense matter equation of state. These observations eliminated 
whole classes of models. Combined with the recent observations of massive neu¬ 
tron stars, the attractive nature of the neutron-neutron interaction at low mo¬ 
menta means the equation of state must be soft at low density with a rapid 
transition to a high-pressure when the higher-momentum neutron-neutron and 
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three-neutron interactions become important. 

Astrophysical observations are providing better constraints on the whole mass- 
radius relation of dense matter, which is critical to a real understanding of neutron 


stars (11,12,13,14). In addition, neutron-star mergers are increasingly seen as 


an important site for heavy-element synthesis, providing even more important 
reasons for understanding the properties of neutron-rich matter. In the near 
future, we can look forward to the observation of gravitational waves from these 


neutron-star mergers (15), providing the most direct evidence of the structure of 


neutron matter and neutron stars. 

2 VERY LOW DENSITIES: NEUTRON MATTER AND COLD 
ATOMS 


At very low densities the equation of state of neutron matter is determined by the 
s-wave (spin-0, isospin-1) neutron-neutron interaction. This interaction is very 
attractive, almost large enough to produce a bound di-neutron. The neutron- 
neutron scattering length is -18.5 fm much larger than typical nuclear 

scales; the neutron-neutron effective range is about 2.7 fm. 

Bertsch had proposed a model of low-density neutron matter that was simply a 


zero-range interaction tuned to infinite scattering length (18), now called the uni¬ 


tary limit in cold atom systems. Remarkably, experiments in cold atom systems 
shortly thereafter provided constraints on the properties of the unitary Fermi gas 
(and therefore also indirectly of low-density neutron matter). The relevant pa¬ 
rameter in these strongly interacting low-density systems is the Fermi momentum 
kF = (3-7T 2 p) 1//3 , where p is the number density, times the scattering length a; in 
the limit that the scattering length goes to infinity all properties of the system 
are fundamental constants times the corresponding free Fermi gas quantity. The 
equation of state in this limit is characterized by the Bertsch parameter £: 


E = ^ Efg = f I L k 


F, 


( 1 ) 


where E F q is the energy of a two-component free Fermi Gas at number density 
p. Similarly the superfluid pairing gap is characterized by the ratio 5 of the gap 
in the unitary limit to the Fermi energy of the non-interacting system: 

1 h 2 

A = E(N + l)--[E(N) + E(N + 2)] = 5 E F = 5 — k 2 F . (2) 

2 2m 

Several properties of the unitary Fermi gas have been studied both theoretically 
and experimentally. This is a very strongly correlated system, with a pairing gap 
of the order of the Fermi energy (5 ~ 1), unlike traditional superconductors 
(where 6 ~ 10 -4 ). We concentrate in this section on comparisons of the equation 
of state and superfluid pairing gaps of cold atoms and neutron matter. While 
neutron matter is much more complicated, with a significant effective range and 
spin-dependent interactions that give small but finite corrections at low density, 
important similarities remain between the two systems. A detailed review of 


pairing in neutron matter can be found in (19), while comparisons to cold atom 


physics can be found in (20,21). 
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2.1 Interaction and Equation of State 


The equation of state is in a sense the most fundamental property of neutron mat¬ 
ter. At low densities, as we have described, it can be related to the properties of 
cold atoms, which have a tunable scattering length and an effective range that 
is nearly zero compared to the average interparticle spacing. Neutron matter, in 
contrast, has a fixed scattering length and effective range. They can be compared 
to cold atoms by comparing systems at the same scattering length times Fermi 
wave number, kpa , as well as the same effective range times Fermi wave num¬ 
ber, fci?r e , (i.e. by using dimensionless parameters characterizing the low-density 
system). 

The equation of state for neutron matter and cold atoms have been compared 
in Quantum Monte Carlo calculations (20,22). For cold atoms we assume a 
nearly zero-range interaction tuned to give the same scattering length times Fermi 
momentum as neutron matter. The details of the interaction are not important 
as long as we are at low density. 

For neutron matter various simple interactions have been considered. The 
simplest is a pure s-wave interaction that acts only in relative s-waves, as these 
are the dominant interactions at low density (see the discussion in Section 3.1 
below). The potential is adjusted to have the same scattering length and effective 
range as that inferred from the neutron-neutron interaction (which is, in general, 
active in higher partial waves). We also consider a more realistic interaction 


that contains both s- and p-wave interactions (23) as the latter start to become 
relevant at higher densities. These simple interactions include all the physics 
required for the study of low-density matter, but at higher densities (and hence 
higher Fermi momenta) the full neutron-neutron and three-neutron interactions 
must be considered. 

At low density neutron matter is a very strongly paired superfluid, requiring 
the inclusion of pairing into a realistic treatment of even the equation of state. 
These calculations calculate the ground state through Quantum Monte Carlo 
methods, specifically Diffusion Monte Carlo. They project out the ground state 
through: 

N 

T 

AC 


I'hn 


= exp = JJexpf-tf —j |T t ), 


2—1 


(3) 


where we take the limit of large r to reach the ground state and is an ap¬ 
proximate trial wave function for the system. In this case we take the trial wave 
function to be of a BCS-Jastrow form: 


|\Fr) 


n/M 

A 

n^) 

}<3 


. *>■? 


(4) 


where <j){r ) is an s-wave pairing function describing the pairing of spin up and spin 
down neutrons, and the Jastrow function / can be used to adjust the short-range 
behavior of the wave function independently of the superfluid pairing. For the 
appropriate choice of 4>(r), this wave function reduces to the standard Jastrow- 
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Figure 1: The equation of state of low-density neutron matter compared to that 
of cold atoms at the same value of Fermi momentum times scattering length 
(fcjra). The left side compares cold atoms and neutron matter (see text), and the 
right panel shows neutron matter results for different methods over a wider range 
of kpa. Figures taken from (21,22). 


Fermi Gas wave function. The radial form of the function <p{r) and f(r) are 
determined in variational calculations. 

These calculations have a fixed-node approximation that implies they provide 
variational upper bounds to the true energy. They have proven to be very accuate 
in studies of cold atom systems, where accurate lattice calculations without a 
fixed-node approximation are available (©. These calculations are also in very 
good agreement with cold atom experiments |2). 

The results of the equation-of-state calculations are shown in Fig. [l] The 
left panel compares neutron matters and cold atoms at very low density. The 
vertical axis indicates the ratio of fully interacting energy to the energy of the free 
Fermi gas at the same density, the horizontal axis is the Fermi momentum times 
the scattering length kpa ; on the upper axis the equivalent Fermi momentum 
for neutron matter is indicated. At extremely low densities, or equivalently small 
value of kpa, analytic results are available ( 24|25 ), and the higher-order Lee-Yang 
result is plotted as a line in the figure. 

Results for cold atoms with zero effective range are plotted as filled blue circles, 
in the limit of infinite kpa these should approach 0.37. Cold atom results for the 
dependence on the effective range are also available, the equation of state can be 
expanded in terms of kp: 


E /E fg = £ + S kp r e + ..., 


(5) 


where S = 0.12(3) is a universal constant that has been determined in the lattice 
calculations and in Diffusion Monte Carlo (3,26). Using the above equation of 
including the experimental neutron-neutron effective range r e gives the dashed 
line in the figure. 
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Figure 2: The pairing gap of low-density neutron matter compared to that of 
cold atoms at the same value of Fermi momentum times scattering length (left). 
Comparison of different calculations of the neutron matter gap (right). Figures 
taken from (20,22). 


The calculations of neutron matter with the spin singlet s-wave interaction 
from AV18 gives the solid red points, correcting the p-wave interactions from 
AV4 gives the green squares. At very low densities all these calculations are very 
similar. At slightly higher densities the correction from the p-wave interaction 
is slightly repulsive, and is in good agreement with the effective range expansion 
above. 

The right panel shows a comparison of various methods for the neutron matter 
equation of state over a somewhat wider range of densities. Methods include 


Fermi Hypernetted chain resummation techniques (27,28), the Brueckner-Bethe 


Goldstone expansion (29), effective field theory (30) and several Quantum Monte 


Carlo methods including Green’s function Monte Carlo (31) and Auxiliary field 


diffusion Monte Carlo (32,33). All the calculations are in reasonable agreement, 


indicating a soft neutron matter equation of state at low density. 


2.2 Superfluid Pairing 

Neutron pairing at low density is important in both neutron-rich nuclei and in the 
crust of neutron stars. Pairing in nuclei and matter has been a long-studied topic, 


see a review by Dean and Hjorth-Jensen (34). For pairing at low density, recent 


work in cold atom systems, both theoretical and experimental, has advanced our 
understanding of pairing in the strongly superfluid regime. Experiments and 
calculations indicate that the pairing gap in neutron matter is quite substantial, 
reaching a peak of approximately 30 per cent of the Fermi energy. This is the 
largest Fermion pairing gap known in nature, and only slightly smaller than the 
45 per cent pairing found in cold atoms at unitarity (j4j[5„6). 

The pairing gap for matter can be calculated by computing the energy of 
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even and odd particle systems and extrapolating to the thermodynamic limit. 
For traditional superfluids and superconductors this is extremely difficult as the 
small value of the pairing gap implies a very large coherence length, in principle 
requiring extremely large simulations with very small differences in energy. For 
this reason Quantum Monte Carlo methods had not previously been used to 
extract superfluid pairing gaps. 

Analytic techniques, are, of course available. In nuclear physics the BCS equa¬ 
tions have often been used to study pairing in neutron matter. At low densities 
the s-wave pairing predictions from BCS theory are mostly independent of nu¬ 
clear interaction models because the interaction is very well constrained at low 
momenta. However corrections to this mean-field theory varied widely, as dis¬ 
cussed in (34). 

The pairing gap of cold atoms and neutron matter are compared in the left 
panel of Fig. [2j At very low kpa the Gorkov-Melik-Barkhudarov analytic correc¬ 
tions to BCS theory are available (35). These corrections reduce the pairing gap 
by a factor of l/(4e) 1 / 3 from that obtained in the simple BCS theory, resulting 
in a pairing gap of: 


A °(M = 


8 


(4 e ) 1 /3 e 2 2m 


exp 


(—) 

\ 2 akF ) 


( 6 ) 


This pairing gap is, however, accurate only at small values of fcpa, and hence not 
applicable near the peak of the pairing gap. 

At these larger values it is again useful to compare with the BCS-BEC transi¬ 
tion in cold atom systems. In the extreme BEC limit, the BCS equation yields 
a pairing gap of 1/2 the binding energy of a pair, as the unpaired fermion can 
propogate nearly freely in the medium when the pairs are very strongly bound. 
There will be no further polarization corrections to this limit. It is natural to 
assume a smooth transition between the BCS limit, where the beyond-mean- 
field limits reduce the gap significantly, to the strongly-paired limit where the 
corrections are very small. 

This smooth transition is in fact observed in cold atom experiments and in 
calculations of cold atoms and neutron matter. The pairing gap has been studied 
extensively experimentally, using both the density harmonically trapped up and 
down spins in a polarized system(6) and by studying the radio-frequency (RF) 
response (36). Analyses of the experimental result indicate a pairing gap of 
0.45 ± 0.05 Ep, or approximately 65 per cent of the pairing gap found in the BCS 
theory. 

The calculated pairing gaps in neutron matter are not as large as those in 
cold atoms because of the finite effective range. Eventually at larger momenta 
the s-wave interaction becomes repulsive. These higher order corrections result 
in a peak of the gap of about 0.3 times the Fermi energy at kpa ~ —5, or 
k F K, 0.25-0.3 fm" 1 . 

A variety of calculations of the pairing gap are shown in the right panel of 
Fig. [2j The simple BCS theory is shown as a solid line, other calculations in¬ 
clude correlated-basis-function calculations (37,38), polarization-potential (39), 
medium-polarization calculation (40), Brueckner (41), renormalization group (42), 
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determinantal lattice Monte Carlo calculations (43), Diffusion Monte Carlo (22), 
and AFDMC results (32,44). 

The pairing gaps do not agree nearly as well as the equation of state. However 
most of the approaches employed are not designed to handle the strong pairing 
that is present for neutron matter. Only the QMC approaches predicted the 
pairing gap in cold atoms and smoothly interpolate between BCS and BEC limits. 
They all result in rather large pairing gaps, with a modest reduction of the BCS 
result for pairing. The behavior of the pairing gap at higher densities is certainly 
interesting and important, particularly the p-wave gap, but requires additional 
study for a reliable prediction. 


3 MODERATE DENSITIES: INTERACTIONS AND EQUA¬ 
TION OF STATE 

Low-density neutron matter, as seen in the previous section, can be described 
by a neutron-neutron potential that is central, i.e. a function of the interparticle 
spacing alone, V(r). As we saw, the essential features of the interaction can be 
captured by the s-wave scattering length a and the effective range r e , i.e. only 
two numbers. 

As the density increases different partial waves and spin-states come into the 
picture. That means that any formulation of the nucleon-nucleon potential needs 
to account for all of the allowed ways in which nucleons can interact with each 
other. Similarly, any quantum many-body method that is used to address in¬ 
teracting nucleons had better be able to handle the different partial waves and 
spin-states (and for the general case also isospin-states). In this section we first 
go over the basics of the operator structure of standard nuclear forces (both phe¬ 
nomenological ones and those following from chiral Effective Field Theory) before 
discussing results for the pure neutron matter equation-of-state (energy versus 
density) resulting from a variety of many-body methods. 

3.1 Neutron-Neutron Interaction 

Nucleon-nucleon (to be precise, np) scattering in different channels can be ex¬ 
perimentally probed, so any nuclear force needs to reproduce a large number 
of scattering phase shifts. The different “channels” in which phase shifts are 
extracted from experiment are conveniently encapsulated in the spectroscopic 
notation 2S+1 Lj : L denotes the orbital angular momentum using S,P, D, and so 
on. The 25 + 1 exponent for two neutrons can be either 2-0 + 1 = 1 (spin-singlet) 
or 2 • 1 + 1 = 3 (spin-triplet). Thus, the low-density s-wave scattering mentioned 
above which is dominant at low densities in neutron matter is given as 1 So- This 
is clearly seen in Fig. [3j which displays scattering phase shifts vs nucleon wave 
number/momentum: at vanishing momenta (corresponding to vanishing densi¬ 
ties) the only phase shift that is non-zero is precisely the one in the 1 So channel, 
(the 3 Si channel, familiar from the deuteron, does not come into play for the 
neutron-neutron case). We remind the reader that a positive phase shift corre¬ 
sponds to an attractive interaction: this implies that at very low densities the 
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Figure 3: np scattering phase shifts vs the wave number of a nucleon in the 
center-of-mass frame (based on the Nijmegen93 partial-wave analysis). Charge- 
independence and charge-symmetry breaking terms would slighly change these 
results for the case of neutron-neutron scattering. Figure taken from (19). 


interaction in the 1 So channel is attractive, turning repulsive at around k = 1.7 
fm” 1 (thereby retroactively explaining the gap closure we saw in the previous sec¬ 
tion). At higher densities/momenta, 3 P 2 is clearly the dominant channel, leading 
to 3 P 2 _3 F 2 pairing (the coupling being due to the tensor operator, see below). 

Traditionally, the effective reproduction of phase shifts was accomplished via 
the breakup of the nucleon-nucleon potential in different channels, exemplified 


by the Argonne family of potentials (45). Here is the Argonne v8’ (AV8’) (23) 


case: 


Fi2 = 5> p (r)0? 2 , 

p =i 

Here 0\ 2 are spin-isospin dependent operators, which are given by: 

Oi 2 ~ 1,8 = (1, cr i • cr 2 , S 12 , L • S) x (1, n • r 2 ) ■ 


(7) 


( 8 ) 


In this expression Si 2 = 3cri ■ r (72 • r — cti • (T 2 is the tensor operator, where 
r = ri — r 2 is the nucleon separation vector. Similarly, L = — ir x (Vi — V 2 ) is 
the relative angular momentum and S = cri + <T 2 is the total spin for the pair 
12. An nn pair can only exist in an isotriplet state (T = 1), so t\ • = 1. 

The v p (r) functions in Eq. Q contain one-pion exchange at large distances and 
high-quality phenomenology at intermediate and short distances. 

More recently, the nuclear physics communicty has embraced chiral Effective 
Field Theory (EFT) interactions, in an attempt to connect with the symmetries 
of the underlying fundamental theory (Quantum Chromodynamics) (46,47). In 


what follows, we provide a bare-bones description of chiral EFT, starting with 
the features that are common to all modern chiral EFT potentials. Chiral EFT 
interactions employ a separation of scales, that between pion and vector meson 
masses, and attempt to systematically expand in a small parameter that is the 
ratio between the two. Assuming the power counting employed is self-consistent, 
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this provides a hierarchy of forces controlled by the power of the expansion. Thus, 
chiral EFT systematically includes the analytically derived long-range one-pion 
exchange, as well as the intermediate-range two-pion exchanges, and so on. In 
addition to this, the power counting naturally leads to consistent three-nucleon 
and many-nucleon forces (e.g. at leading-order, LO, and next-to-leading order, 
NLO, no 3N forces are predicted). The soft scale, identified with the pion mass 
above, is generally taken to be that of a relevant momentum scale. In addition to 
the known pion exchanges, chiral EFT interactions also include phenomenolog¬ 
ical short-range terms, typically written down as Dirac-delta functions/contact 
terms: these do not reflect any deeper physics other than respecting all relevant 
symmetries and the power counting. 

Before we look at specific expressions for the different parts of a chiral EFT 
interaction, we go over some definitions. Since chiral EFT interactions are based 
on a low-momentum expansion, it was natural that they were first formulated 
in momentum space. If we denote the incoming and outgoing relative momenta 
by p = (pi — p 2 )/2 and p 7 = (p) — p 2 )/2, respectively, then we can define the 
momentum transfer q = p 7 — p and the momentum transfer in the exchange 
channel k = (p 7 + p)/2. Generally speaking, terms that depend on q are local 
and those depending on k are non-local. The power expansion mentioned above 
is therefore an expansion in powers of q and k. Since our discussion is meant to 
be somewhat pedagogical, we limit ourselves below to the first three orders in the 
chiral expansion though, obviously, the pion exchanges and contact terms have 
also been extensively studied at higher orders. 

Let us now look at the one- and two-pion exchanges in momentum space: 

V™Eo = u t(q 0 (ri • r 2 ) Oi • q) (<r 2 • q) 

V ™NLO = W C (q) (ri • r 2 ) + Vs(q) (<n ■ <r 2 ) + V T (q) (oq ■ q) (<r 2 • q) 

V ™n 2 LO = V c(q) + Ws(q) (n • r 2 ) (oq • <r 2 ) + W T {q) (ri • r 2 ) (oq • q) (oq • q) 

(9) 

Note that all three of these expressions are functions of q alone (i.e. they are 
not functions of k). Thus, they can be trivially Fourier transformed and take the 
following form in coordinate space: 

Us(r) (ri • r 2 ) (<ri • oq) + U T (r) (ri • r 2 ) S i2 
Wc(r) (ri • t 2 ) + Vs(r) (oq • cr 2 ) + V T (r) S 12 
V c (r) + W s (r) (ri • r 2 ) (oq • cr 2 ) + W T (r) (n • r 2 ) S 12 (10) 

where we are only showing the finite-range parts of the pion exchange terms. It 
is easy to see, for example by focusing on the one-pion exchange term, that the 
(<Ti • q) (<t 2 • q) term in momentum space gives rise to cr\ • er 2 and ,SV 2 terms in 
coordinate space. 

Note that the above expressions for the pion exchanges do not contain relativis¬ 
tic 1 /m 2 N and 1/mjv corrections for the one- and two-pion exchange, respectively 
(where mjv is the nucleon mass). These would also involve (isoscalar and isovec¬ 
tor) spin-orbit terms. This means that in chiral EFT the only spin-orbit term up 


T/coord _ 

V l7 t,LO ~ 

t/ coord 
v 2n,NLO ~ 

T/coord _ 

V 2tt,N 2 LO ~ 
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to order N 2 LO appears in the short-range contacts, which, as mentioned above, 
respect all relevant symmetries as well as the power counting. Qualitatively, 
working in momentum space, this means that at LO all 4 possible spin-isospin 
combinations should be present, without any momentum dependence. At NLO, 
something analogous holds: there are 14 different contact interactions that all 
contain a momentum-squared in some form. This can be summarized as: 


K™n™LO = Q i + a 2 tTi ■ cr 2 + a 3 ri • t 2 + a 4 <ri • <r 2 n • r 2 
V mngNLO =71 <T + 72 q 2 CTi ■ <T 2 + 73 Q 2 D ■ T 2 + 74 tfVl • <T 2 Ti • T 2 

+ 75 k 2 + 76 A ; 2 07 • <r 2 T 77 A; 2 ri • t 2 + 73 A: 2 ^! • cr 2 n • t 2 

+ 79 (o"i + cr 2 )(q x k) + 710 (07 + cr 2 )(q x k)n • r 2 
+ 7n( cr i • q)(^2 • q) + 712(0-1 • q)(o- 2 • q)n • r 2 

+ 713(0-1 • k)(cr 2 ■ k) + 714(0-1 • k)(cr 2 • k)ri • r 2 (11) 


(This is in momentum space, leading to a Dirac delta 5(r) in coordinate space, 
hence the name contacts.) It’s obvious that this general set of contacts contains 
both q and k on an equal footing. In reality, only 2 and 7 of these contacts at LO 
and NLO, respectively, are linearly independent. As a rough mapping between 
momentum space and coordinate space, we note that the isoscalar spin-orbit term 
(cr 1 + cr 2 )(q x k) corresponds to L • S. 

All the chiral EFT terms mentioned above need to somehow be cut off (“reg¬ 
ulated”) at high momenta in order to avoid infinities. Thus, they are multiplied 
with a regulator function: 


fk{p',p) = exp 



( 12 ) 


where n is some power and A is known as the cutoff parameter, typically taken to 
be 400-500 MeV. This regulator and cutoff are used in the standard momentum- 
space formulations of chiral EFT by Entem-Machleidt (EM) (48) and Epelbaum- 
Glockle-Meissner (EGM) (49). In addition to this, the EGM approach also uses a 
spectral-function regularization on the two-pion exchange terms given above (49), 
leading to a 2nd cutoff parameter, A which is typically taken to be larger than 
A (700 MeV and above). It is also significant that nuclear forces with explicit A 
degrees of freedom are being investigated. (50,51) 


More recently, two new chiral EFT potentials have appeared: i) N 2 LO op t (52), 
which carried out a high-quality re-fit of the low-energy constants appearing in 
the contact terms, like those shown above, as well as the pion-nucleon couplings, 
and ii) local chiral EFT, (53|54) which made use of the afore-nrentioned ambiguity 
in selecting the contacts such that they are quasi-local to begin with. [The latter 
approach also made use of a local regulator (i.e. one that is a function of q or r), 
as also adopted by a very recent mixed version of chiral EFT (55), which employs 
a local regulator and non-local contacts.] As will be seen in the many-neutron 
results discussed in the following subsection, the N 2 LO op t potential is very soft, a 
feature which is useful to many-body methods that work in momentum space. On 
the other hand, the local chiral EFT approach has varied the cutoff A, producing 
a family of potentials, ranging from hard to quite soft. 
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3.2 Equation of State 


While the nuclear force employed is an important component in describing neu¬ 
tron matter, just as important is the quantum many-body theory that is used 
to study the many-neutron problem. Details on specific results are given below, 
but let us first go over some general features of nuclear many-body methods. We 
limit our discussion to those methods that have been recently reformulated for 
(or applied to) neutron matter. 

The first large distinction we could make is between ab initio methods, which 
try to solve the many-neutron Schrodinger equation employing varying degrees 
and types of approximations, and phenomenological approaches, which focus on 
reproducing and predicting experimentally relevant properties without necessar¬ 


ily seeking a connection to nuclear forces (56). Ab initio methods, in turn, can 


be grouped according to another general distinction, that between perturbative 
and non-perturbative methods. The prototypical perturbative method employed 
to describe neutron matter is Many-Body Perturbation Theory (MBPT), an ap¬ 
proach which attempts to sum all classes of diagrams up to a specified order in 


interaction lines (57,58,59). On the other hand, non-perturbative nuclear meth¬ 
ods can also be distinguished between a) resummation schemes, which sum a 
specific class of diagrams to all orders in the interaction, typically employing a 


ladder approximation, (60,61), b) Coupled Cluster theory, which works by gener¬ 
ating np — nh excitations of a reference state, typically truncated at the doubles 


level (62), and c) Quantum Monte Carlo methods, in which the many-body prob¬ 


lem is solved stochastically, (63,33,53,64 65). One could legitimately add to this 


list d) density functional theory, which de facto includes complicated many-body 
correlations and is applicable throughout the nuclear chart. 

Given the significance of energy-density functional approaches of heavy nu¬ 
clei, it is worthwhile to examine the interplay between ab initio calculations 
for neutron matter and EDF approaches. Ref. (66) explicitly showed that the 


equation-of-state (EOS) of neutron matter can be used to eliminate as unphysi¬ 
cal several Skyrme energy-density functionals. As a result, the EOS of neutron 
matter has been, for a while now, used as a constraint in energy-density func¬ 
tional theories of nuclei (56). Recently, the interplay of Skyrme functionals and 


pure neutron systems was exploited in the context of the polaron/impurity prob¬ 


lem (67), and more generally (68) - see also section 4 of the present review. 


The idea behind energy-density functional theories of nuclei is to use as input as 
much known physics as possible, in order to make predictions in regions where 
no experiment is possible (yet) and fully microscopic calculations are unrealistic. 
The use of microscopic results for pure neutron properties, then, ties in with 
that agenda, extending the input constraints beyond experiment: it stands to 
reason that neutron-rich nuclei predictions will depend critically on pure neutron 
constraints. Of course, the resulting EDFs can be applied to the context of neu¬ 


tron stars as well (69,70), again addressing systems that are beyond the reach 


of ab initio approaches. Such density functionals could also explicitly include a 
neutron-star crust equation-of-state as input. 

We now go over selected results on the equation of state of pure infinite neutron 
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(a) (b) 


Figure 4: Equation-of-state of intermediate density neutron matter using only 
NN interactions as input. Shown are results for different chiral EFT potentials, 
different orders in the chiral expansion, and different quantum many-body meth¬ 
ods. (a) Quantum Monte Carlo results, (b) All other many-body results. Details 
are provided in the main text. 


matter, following from several quantum many-body approaches. All of these 
results use as input chiral EFT neutron-neutron interactions. It’s important to 
note that if one is limited to the NN sector alone, the chiral EFT approach 
has no built-in advantage over phenomenological approaches like the Argonne 
family of potentials or CD-Bonn. One of the main distinguishing features of 
the chiral EFT approach is that it provides guidance (through its order-by-order 
expansion) on how to consistently pick the form of the 3N interaction. Thus, 
the results summarized in this subsection should be seen as the foundation of an 
approach that combines NN+3N forces, which become even more important at 
higher densities, as discussed later in this review. 

We show in Fig. [4] results for the equation of state of neutron matter, fol¬ 
lowing from calculations that use only nucleon-nucleon interactions as input. 
Specifically, we show MBPT results using local N 2 LO interactions as input (54) 
(specifically the soft 400 MeV potential), MBPT results using an N 3 LO EM 
potential (59), MBPT results using N 2 LO op t (71), self-consistent Green’s func¬ 
tions results using N 2 LO op t (60), coupled-cluster doubles (CCD) results using 
N 2 LO op t (62), AFDMC results using local N 2 LO interactions (54), Auxiliary- 
Field Quantum Monte Carlo (AFQMC) results using an N 3 LO EM potential (65), 
and configuration-interaction Monte Carlo (CIMC) results using N 2 LO op t (64). 

We see that at intermediate densities essentially all many-body methods with 
different input potentials give basically the same answer qualitatively (though 
that ceases to be the case closer to saturation density). At a more detailed level, 
within MBPT Fig. [ 4)3 seems to imply that as we go from N 2 LO to N 3 LO the 
energy decreases. Also, the MBPT plus soft local N 2 LO results are essentially 
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identical to the AFDMC plus soft local N 2 LO ones (54). Intriguingly, the 414 


MeV (soft) EM potential plus AFQMC of Ref. (65) leads to results that are 
similar to the 500 MeV (hard) local chiral potential plus AFDMC of Ref. 


The MBPT results using N 2 LO op t of Ref. (71) underline how soft/perturbative 
the N 2 LO op t potential is: this implies that non-perturbative methods are not 
necessarily required to handle it. A distinguishing feature of the AFDMC results 


of (54) is the use of a fully non-pertubative QMC method to probe both soft 
and hard chiral potentials. Overall, it is clear from Fig. [4] that issues like the 
choice of different regulators and different cutoffs, as well as the non-locality vs 
locality of the potential employed, should be further investigated in detail. At 
vanishingly small densities the effects of pairing (discussed in an earlier section) 
would probably alter this picture, but such differences would likely not be visible 
at this scale. 


3.3 Three Nucleon Interactions 

It is well accepted that modern nucleon-nucleon interactions cannot describe the 
binding energies of nuclei with A>3, and they have to be combined with three- 
body forces. In the framework of chiral EFT, as discussed in Sec. 3.1 the ap¬ 


pearance of three- and higher many-body forces naturally emerges in the chiral 
expansion. However, in other approaches like for the Argonne interactions, the 
contribution of four-body forces is expected to be negligible. This observation 
comes from the fact that in nuclei and matter the three-body potential energy is 
(in magnitude) few percent than the two-body potential energy. The Urbana-IX 
(UIX) potential has been introduced to give a correct description of few-body 
systems. The UIX includes the Fujita-Miyazawa term that describes the p-wave 
exchange of two pions between three nucleons, where the intermediate state has 
one nucleon excited to a A, that is eventually integrated out. This term is the 
longest-range three-nucleon interaction and a very similar term also arises as the 


leading three-nucleon contribution (at N 2 LO) in chiral perturbation theory (46). 


The Urbana models also include a spin-isospin independent short-distance term, 


that models the exchange of more pions between nucleons (72). The UIX force 


was originally proposed in combination with the Argonne AV18 and AV8' (73). 


It has been constructed to correctly describe the saturation of nuclear matter 
at po = 0.16 fm -3 , and to reproduce the binding energy of the 4 He. Although 
this model of the three-nucleon force does not fully alleviate the underbinding 
problem in light nuclei, it has been extensively used to study the equation of 


state of nuclear and neutron matter (74,33,75,76). It has to be noticed that the 


AV18+UIX Hamiltonian supported heavy mass neutron stars years before they 
have been observed. 

Other phenomenological models of three-body force, such as the Illinois forces, 


have been obtained by fitting the binding energies of light nuclei up to A<8 (77), 


and they give an excellent description of both ground- and excited-states of nuclei 
up to 12 C (78). In addition to the Fujita-Miyazawa term the Illinois models also 
include s-wave two-pion-exchange and also three-pion-exchange ring diagrams. 
Unfortunately, the three-pion rings included in the Illinois forces leads to a very 
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strong attraction in pure neutron-matter (79,80) at high densities, and produces 


EOS too soft to be compatible with observed neutron star masses and radii. This 
is a sign of important deficiencies in the model of the three-nucleon forces, which 
will become especially important at higher density and larger isospin asymmetry. 
Hamiltonians based on chiral EFT would be very valuable to address these issues. 
However, such interactions should be constructed to reproduce nucleon-nucleon 
scattering data up to large energies. A qualitative estimate is given by observing 
that the momentum of two neutrons is given by k ~ y/Elab nr/2, and then if 
this is compared to the Fermi momentum of homogeneous matter we have p ~ 
(Eiab m/2) 3 / 2 /3 it 2 (of course, this is only a back-of-the-envelope calculation, so 
it does not account for possible prefactors). This clearly shows that if one wants 
to describe important correlations in neutron matter up to large densities, then 
a nucleon-nucleon interaction that accurately fitting phase-shifts up to large 
will be important. Of course, this raises the questions of which functional form 
for the three-nucleon force should be employed, as well as which observables 
the three-nucleon parameters should be fit to, when one interested in describing 
matter at high density. 

4 INHOMOGENEOUS NEUTRON MATTER 


The neutron matter in the exterior of neutron-rich nuclei or the inner crust of 
neutron stars is not homogeneous, so it is important to study the properties of 
neutron matter in confined geometries and with density gradients. Typically, 
nuclear density functionals have isoscalar and isovector gradient terms that al¬ 
low for a different energy cost of density gradients between isospin symmetric 
nuclear matter and pure neutron matter. The isoscalar gradients can be very 
tightly constrained by the properties of finite nuclei, but the gradients in pure 
neutron matter are much less tightly constrained. To some degree this is a result 
of the fact that we have direct experimental extractions of the charge density 
through electron scattering, but less direct knowledge of the neutron density. At 
present there is a large extrapolation of experimental knowledge to the nearly 
pure neutron matter found between nuclei in the inner crust of neutron stars. 

One avenue to extract additional information is to use microscopic theories 
of neutron matter in an external field. Different external fields can be used to 
probe different properties of neutron matter; to date finite systems have been 
studied using harmonic oscillator or Woods-Saxon external one-body potentials. 
These serve to confine the system and to constrain the densities. Originally these 


calculations were limited to very small systems (81), but more recently much 


larger systems have been treated where comparisons to density functionals are 


more meaningful (82). 


Other properties of density functionals can also be studied in the pure neutron 
limit. The relative importance of pairing and shell gaps is tightly constrained by 
the measured properties of nuclei. Shell gaps in neutron matter are expected to 
be much smaller because of the strong pairing. Spin-orbit splitting can also be 
examined by studying systems away from closed shells. 
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4.1 Neutron Matter and Density Functionals 

It is again instructive to compare results for confined systems in cold atoms with 
those in neutron matter. Cold atoms at unitary have essentially no closed shells, 
as discussed below, and hence can be described as an expansion around the local 
density: 


£2 = J v ext {r)p(r) +^(37T 2 ) 2/3 p 5/3 +c 2 Vp 1/2 • Vp 1/2 + c 4 


vy/ 2 vy / 2 


p 


2/3 


+ ..., 

(13) 

where the first two terms provide the simple local-density approximation (LDA). 

The form of this energy density functional is severely restricted by the scale- 
invariance of the cold atom system. All terms in the energy density must scale 
like p 5 / 3 . Such a functional will not work for cold atoms at small kpa as in this 
BCS regime individual particles can propagate over the entire system and the 
fermionic closed shells are very important. While it is not obvious a priori that 
this simple density functional will work at unitarity, microscopic results using 
both Diffusion Monte Carlo and lattice Auxiliary Field Monte Carlo justify this 


simple local density plus gradient expansion form (83). 


Clearly such a simple form will not work well for nuclei or neutron drops, where 
the fermionic shell closures play a critical role. Of course Skyrme and Hartree- 
Fock-Bogoliubov models play a critical role in the study of medium and heavy 
nuclei, including particularly neutron-rich nuclei (8). These density functionals 
have parameters that are fit to the properties of nuclei, and sometimes to the 
calculated bulk properties of neutron matter. A similar strategy should be helpful 
in constraining the gradient, spin-orbit, and pairing terms in very neutron-rich 
nuclei. 


4.2 Neutron Drops 

Here we summarize results for neutron drops in Hamonic Oscillator and Wood- 
Saxon external fields, and compare them to available results for cold atoms. A 
variety of different interactions and methods have been used for calculations of 
neutron drops ( |82[[84 ,.85||86). In Fig. [5] we show results for neutrons trapped 
in a harmonic oscillator trap with a frequency of 10 MeV. Results for lower and 
higher trap frequencies have also been reported. For cold atoms lattice methods 
are used which are free of a sign problem; note that the cold atom results show 
essentially no evidence of the closed shells that would appear for non-interacting 
fermions with N=8, 20, 40, etc. The local density approximation with the bulk 
value of £ = 0.37 for cold atoms is shown as a solid black line, finite drops rapidly 
and smoothly approach this bulk limit. Gradient corrections have been studied in 


Ref. (83). The dashed line indicates the local density approximation for £ = 0.50, 


a value more appropriate for the neutron equation of state. 

Calculations from 8 to 50 neutrons using the AV8’ NN interaction and the AV8’ 
+ UIX NNN interaction are shown. The results plotted here use the AFDMC 
method; results for N up to 16 with GFMC are essentially identical. The three- 
nucleon interaction is not very important for small N, as the density is relatively 
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Figure 5: Energies for neutrons trapped in a harmonic oscillator with uj = 10 
MeV for different interactions and methods compared to Monte Carlo calculations 
for cold atoms at unitarity (open symbols) and two local density approximations 
(solid and dashed lines). In the right panel results for the AV8’ + UIX interaction 
are compared to several pre-existing density functionals. 


modest. For larger N, though, the central density in a harmonic trap can be 
quite large and the three-nucleon interaction makes a considerable difference. 
Results are also shown for no-core shell model (NCSM) calculations using chiral 
interactions and the JISP16 NN interaction. The chiral interactions are very 
similar to the AV8’ results without a three-nucleon interaction; the energies for 
all interactions and methods agree within a few percent. 

These calculations show significant shell closures, in contrast to the cold atom 
results. These shell closures arise because of the effective range in the neutron- 
neutron interaction, and the concominant smaller size of the superfluid pairing 
gap. In cold atoms only superfluid pairs propagate across the whole system, while 
for neutrons the shell closures indicate the single-particle picture still survives to 
some degree. The shell closures are not as strong as in nuclei, however. Pair¬ 
ing gaps are also evident from the odd-even staggering, spin-orbit, and gradient 
corrections have been examined in many of these studies (84,85). 

In the right-hand panel, results for the AV8’+UIX interaction are compared 
with a variety of previously existing density functionals. These funtionals do not 
fully describe the neutron drop results, in general their energy is too low and 
spin-orbit and pairing are not completely correct. Modest modifications to the 
isovector (full neutron) gradients, pairing and spin orbit give density functional 
results in the full black line, labeled SLY4-adj. These do a good job of repro¬ 
ducing results in both 5 and 10 MeV frequency traps except for the smallest 
systems considered. More recently it has been shown that new density function¬ 
als, i.e. UNEDF0 (87) and UNEDF1 (88), can be created that simultaneously 
fit nuclei with accuracy comparable or better than the existing functionals, and 
also reproduce the neutron drop results. 

In figure [6] we show results in a Wood-Saxon well. These more nearly mimic 
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(b) 


Figure 6: Energies for neutrons trapped in a Wood-Saxon external field compared 
with density functionals (left panel) and radii of neutron drops in HO and WS 
wells compared to density functionals (right panel). 


nuclei since the central density is rather fixed, and hence they are not so sensitive 
to the high-density behavior of the equation of state. The same modified density 
functional that fits the harmonic oscillator results also reproduces the Wood- 
Saxon results, and simultaneously describes the calculated radii of drops in both 
external fields. 

Further studies of inhomogeneous neutron matter are warranted, including 
effects of the three-nucleon interaction and different external fields. In the inner 
crust of a neutron star, for example, the neutrons are not small finite systems. 
Therefore it would be very useful to constrain the gradient terms with periodic 
external fields with a range of momenta, and see if the same density functional 
can describe these systems as well as finite drops. 

5 HIGHER DENSITIES: NEUTRON MATTER AND NEUTRON 
STARS 

At higher densities the three-nucleon interaction becomes even more important. 
These play an important role in the mass/radius relation for neutron stars. In this 
section we describe the equations of state obtained at higher densities including 
three-nucleon interaction results and the equation of state, describe the relation 
between the symmetry energy E sym , the difference between neutron matter energy 
and isospin-symmetric nuclear matter energy at nuclear saturation density, and 
its density derivative L. Finally we describe the impact of these constraints on 
the neutron star mass-radius relation as well as describe briefly the importance 
of three-hadron neutron-neutron-lambda interactions in the equation of state of 
high-density neutron matter. 

5.1 Equation of State 

The EOS of neutron matter around nuclear densities is important for several 
reasons. As outlined in the previous sections, in this regime several channels 
contributing to the nucleon-nucleon interaction are important, and thus the EOS 
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Figure 7: Equation-of-state of intermediate density neutron matter using both 
NN interactions and 3NF as input. Shown are results for different chiral EFT 
potentials, different orders in the chiral expansion, and different quantum rnany- 
body methods, (a) The standard APR results, along with Coupled Cluster and 
Auxiliary-Field Quantum Monte Carlo results (b) Many-body perturbation the¬ 
ory, Self-Consistent Green’s Function, and pp ladder approximation results. De¬ 
tails on the interactions employed are given in the main text. 


provides a direct comparison between different frameworks for the nuclear Hamil¬ 
tonian. Although the maximum mass of neutron stars is dominated by the EOS 
at very high densities, their radius is determined by the pressure in the region of 
about 1 to 2 po) and then measurements of radii of neutron stars can be used to 
constrain the EOS. 

We show in Fig. [7] results for the equation of state of neutron matter at some¬ 
what smaller densities, following from calculations that use both nucleon-nucleon 
and three-nucleon interactions as input. Specifically, we show we show MBPT re¬ 
sults using EM and EGM N 2 LO interactions as input (the lower end of this band 
is shown using a dashed line), (58) MBPT results using N 3 LO EM potentials, 
(58) results following from a particle-particle ladder approximation using N 2 LO 
EM potentials, (61) as well as pp ladder results using an N 2 LO EM potential, 
self-consistent Green’s functions (SCGF) results using N 2 LO op t, (60) as well 


as SCGF results using EM N 3 LO as input, (60) Coupled-Cluster with doubles 
and including some triples effects, CCD(T), results using the N 2 LO opt interaction, 
(62) Auxiliary-Field Quantum Monte Carlo (AFQMC) results using an N 3 LO EM 
potential, (65) as well as the frequently cited Akmal-Pandharipande-Ravenhall 
results using AV18 plus UIX ( |74| ). 

The most obvious feature of Fig. [7] in comparison to Fig. [4] is that three-nucleon 
interactions in neutron matter are repulsive, i.e. the energy goes up: this is in 
contradistinction with what holds in light and medium-mass nuclei. Another 
general aspect of this figure relates to the transition from N 2 LO to N 3 LO: in 
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Fig- [Zb for MBPT, SCGF and the pp ladder approximation results are shown for 
both orders in the chiral expansion, and in both cases the energy decreases. This 
is despite the fact that the MBPT results use consistent interactions, i.e. include 
an N 3 LO 3N force along with the N 3 LO NN force, while the pp ladder and SCGF 
values, like most other many-body results currently available, combine the N 3 LO 
NN force with N 2 LO at the 3N level. Also prominent is the fact that the SCGF 
N 2 LO op t results reasonably match the EM N 2 LO pp ladder values (and similarly 
for the corresponding EM N 3 LO ones), which is probably not unexpected (with 
the exception of the low-density region). A final noteworthy feature is that the 
CCD(T) results appear to match rather nicely with the APR values. A conclusion 
that the Figs. [4] & [7] really underline is the significance of producing error bands 
using more than one potential as input, if possible employing a non-perturbative 
or otherwise controlled approximation scheme. 

Around density po, the EOS of pure neutron matter is directly related to the 
symmetry energy and its slope, and the analogous (and much more challenging) 
calculation of symmetric nuclear matter is not needed. The symmetry energy 
Esym(p) is generally defined as the difference between the energy of pure neutron 
matter and the energy of symmetric nuclear matter with the same total barion 
density. In terms of the isospin asymmetry, x = p p /p, the energy per nucleon of 
isospin asymmetric nuclear matter is commonly expanded in even powers of x, 

E(p, x) = E 0 (p ) + Egl(p)(l - 2.x) 2 + Egl( 1 - 2x) 4 + ... , (14) 

where E is the energy per nucleon of the system, Eq(p) = E(p,x = 0.5) is the 
EOS of symmetric nuclear matter, and Esym(p) and higher order corrections are 
ignored here. The symmetry energy E s ym is then given by 

E S ym{p) = E(p, 0) - E 0 (p) , (15) 


where E(p, 0) is the EOS of pure neutron matter. Near the nuclear saturation 
density, po, there are a number of constraints on the EOS of infinite nuclear matter 
from nuclear masses, charge radii, and giant resonances. The extrapolation of 
the binding energy of heavy nuclei to the thermodynamic limit yields Eq(po) = 
— 16.0 ±0.1 MeV (89). Because the pressure is zero at saturation, the symmetry 
energy can be expanded around saturation as 


E sym (p) 


— Es ym ± 


PO 


L p- po 
3 Po 


(16) 


where E sym is the symmetry energy at saturation, and L is a parameter related 
to its slope. It is convenient to first focus to E sym at po for several reasons. The 
calculation of symmetric nuclear matter is much more challenging than pure neu¬ 
tron matter, mainly because the Hamiltonian is more complicated, and isospin 
symmetric matter is more strongly correlated through the tensor interaction act¬ 
ing in S- and D- waves. In particular, in the low density regions clustering effects 
are very difficult to be included in the calculation of the EOS, and/or isolated 
from experimental measurements. 
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Using QMC methods, the EOS of pure neutron matter is typically calculated 
with a simulation of 66 particles in a box imposing periodic boundary condi¬ 


tions, for which finite size effects are generally very well under control (79,33). A 


comprehensive study of the model of the three-neutron interaction has been pre¬ 
sented in Ref. (75), with particular emphases to the role played by short range 


correlations. Within the Urbana/Illinois models, the main contribution of the 
three-body force is given by the short-range part, whose structure in the UIX 
model is 

V-jk = ArY^ T 2 (m^r ij)T 2 (m^rjk ), (17) 

eye 

where is the pion mass, and 


T( x ) — ( 1 -|-1 2 

x ar 


P ~ x 

—e(x) 

X 


(18) 


The function £(x) = 1 — e~ cx is a cutoff function to regularize T(x) at small 
distances. The V R term basically models the four-pion exchange between neu¬ 
trons (90). In order to address the role of the above short-distance behavior, 
other forms for V R k have been considered, in addition to different models of 
intermediate- and long-range contributions of the three-body force. These terms, 
where two or three pions are exchanged between neutrons, with the creation of 


A excited states, are described in Ref. (77). 


The resulting EOS are shown in Fig. [8j where we compare the results obtained 
with the AV8 7 and AVS'+UIX Hamiltonians with several EOS obtained using 
different models of three-neutron interactions, adjusted to give the value for E sym 
indicated in the figure. The range of E sym is compatible with several experimental 
measurements (91,92). To obtain a lower symmetry energy with the AV8' model 


of two-body force would require an attractive contribution of the three-body force. 


The three-pion rings give attraction in pure neutron matter (80), but they would 


require a strong short-range repulsion to make the EOS stiff enough to support 
astrophysical observations. 

The EOS calculated using QMC can be conveniently parametrized using the 
following functional form: 


E(p n ) = a ( — 

\Po 


+ M 

HoJ 


ft 


(19) 


where E(p n ) is the energy per neutron as a function of the neutron density p n , 
and the parameters a, a , b, and /3 are obtained by fitting the QMC results. The 
parametrization of the equations of state obtained with the AV8' and AVS'+UIX 


Hamiltonians are reported in (93). 

Using Eq. (16) the value of E sym and L can be easily extracted from the calcu¬ 


lated EOS. The result are shown in Fig. [9], where we compare the results obtained 
using the AV8 7 and AV^+UIX Hamiltonians (red and black symbols), the vari¬ 
ous EOS giving the indicated E sym obtained by changing the three-neutron force 
model (green and blue symbols), and results obtained using the Illinois model of 
three-neutron force that includes three-pion rings where we have independently 
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Figure 8: The equation of state of neutron matter obtained by using various 
models of three-neutron force as described in the text. For each model we impose 
that the energy at saturation is 17.7(1) MeV (blue band), or 16.0(1) MeV (green 
band). The results are compared with the equations of state obtained with the 
AV8' and AVS'+UIX Hamiltonians. In the legend we indicate the corresponding 


symmetry energy at saturation. Figure taken from (93). 
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Figure 9: The value of L as a function of E sym obtained from various EOS. The 
green and blue points with error bars correspond to the various EOS indicated 
by the two colored areas of Fig. [8j and red and black points show the results 
obtained using a two-body force alone and combined with the UIX model. The 
square symbols correspond to results obtained by independently changing the 
cutoff parameters entering in Vr and in the three-pion rings of the three-neutron 


force. Figure taken from (93) 
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changed the cutoff of the intermediate- and short-range part. It is clear that 
within this model the correlation between L and E s ym is very strong. 


5.2 Neutron Star Mass/Radius 

The neutron star matter is mainly composed by neutrons and a few protons. 
When the EOS of the neutron-star matter has been specified, the structure of 
an idealized spherically-symmetric neutron star model can be calculated by inte¬ 
grating the Tolman-Oppenheimer-Volkoff (TOY) equations: 


dP G[m(r) + Airr^P/c 2 ] [e + P/c 2 
dr r[r — 2Gm(r) / c 2 ] 


dm(r) 

dr 


= Auer 1 


( 20 ) 


where m{r) is the gravitational mass enclosed within a radius r, and G is the 
gravitational constant. The above equations are solved by obtaining the energy 
density e and pressure P from the EOS, and by specifying an initial point in the 
integration that is given by the central density of the star. For the EOS with the 


form of Eq. 19 we have 


e = Po 


-rw-) 

Po J \PoJ 
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+ m n | — 
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P = Po 
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( 21 ) 


( 22 ) 


The solution of the TOV equations gives, for a specified central density p c , the 
profiles of p, e and P as functions of radius r, and also the total radius R and 
mass M = m(R). The total radius R is defined to the point where the pressure 
vanishes. The solution of the TOV equations are modified only slightly by mag¬ 
netic fields and temperatures which are expected, and rotation is less than a 10% 
effect for the kinds of M-R curves which we present below. The speed of sound, 
c s in the neutron star interior is e? s = dP/de, and ensuring that this is less than 
the speed of light (and thus the EOS is said to be “causal”) constrains the set 
of possible EOS. Also, the pressure must increase with increasing energy density, 
dP/de > 0, in order to ensure that the neutron star is hydrodynamically stable. 

In practice, the crust of neutron stars cannot be modeled with a pure neutron 
matter EOS. Several approaches to describing the EOS of neutron-star matter 
exist, including the combination of microscopic calculations based on chiral EFT 


with polytropes put forward in (94,95). All the results presented in this section 


have been obtained by using QMC calculations shown in the above sections for 
densities p > p C mst = 0.08 fm -3 . At lower densities the neutron star matter is 


modeled by using the EOS of the crust obtained in Refs. (96,97). When the 


EOS violates causality (this is very often the case of EOS obtained from non- 
relativistic nuclear Hamiltonians), the EOS is switched to the maximally stiff 
EOS. The TOV equations are solved for several values of p c , and the solution 
giving the maximum mass is considered. 















24 


S. Gandotf 



Figure 10: The mass-radius relation for neutron stars based on the QMC neutron 
matter results above. Results are presented for the different EOS given in Fig. [8j 
The numbers indicate the value of E sym of the various EOS. 


By using the EOS obtained from different nuclear Hamiltonian, we can study 
the effect to the neutron star structure. The results of the M-R diagram of 
neutron stars obtained from the EOS calculated in the previous section are shown 
in Fig. [To) Since the radii of neutron stars are almost determined by the EOS 


slightly above po (98), future measurements will provide strong constraints to the 
nuclear Hamiltonian. In particular, radii are directly connected to the pressure 
of neutron matter at po, and then there is a natural correlation between E sym 
and L and radii. In the figure the two bands correspond to the EOS described 
in the previous section (the corresponding values of the symmetry energy are 
also indicated in the figure). The red and black curves correspond to the EOS 
calculated with the AV8 7 two-body interaction alone, and combined with the UIX 
three-neutron potential. The relation between E sym and the radius is evident, 
as the increasing of E sym predicts a neutron star with a larger radius. In the 
figure, the density of the neutron matter inside the star is indicated with the 
orange lines. As anticipated, even at large masses the radius of the neutron star 
is mainly governed by the equation of state of neutron matter between 1 and 2 
po (98). 

As is clear from the figure, the AV8' Hamiltonian alone does not support the 
recent observed neutron star with a mass of 1.97(4)M 0 Q and 2.01(4)M o (10). 
the addition of a three-body force to AV8' can provide sufficient repulsion to be 
consistent with all of the constraints. The results also suggest that the most mod¬ 
ern neutron matter EOS imply a maximum neutron star radius not larger than 
13.5 km, unless a drastic repulsion sets in just above the saturation density (75). 
This rules out EOS with large values of L , typical of Walecka-type mean-field 
models without higher-order meson couplings which can decrease L. We note 
that our analysis suggests it is unlikely that neutron stars have radii lower than 
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10 km (13,14), and/or larger than 15 km, but determining the systematic uncer¬ 
tainties is still an open question (99). 

These recent observations of neutron stars with 1.97(4)Mq (9) and 2.01(4)Mq 
(10) put the most severe constraints on the EOS, although the precise hadronic 
composition is still undetermined. When the EOS is stiff, the corresponding 
chemical potential is large enough that heavier particles, hyperons, can be created 
from neutrons. For example, in the case of neutrons the A particles might form, 
and a fundamental question is about their role to the EOS. Note that other lighter 
hyperons (like S) would form at higher densities because they need that the 
proton fraction is large enough, and consequently their production is suppressed 
compared to A. The main uncertainty in the inclusion of hyperons is because the 
nucleon-hyperon (and the hyperon-hyperon) interaction is not well understood. 
For example, several calculations based on Brueckner-Hartree-Fock suggest that 
when hyperons are included in the EOS, the corresponding mass of neutron stars 
is very low (100). Other calculations based on relativistic mean-field or similar 
methods instead do not exclude that the inclusion of hyperons give neutron stars 
that are not compatible with recent observations. The A-nucleon potential, and in 
particular the importance of the A-neutron-neutron interaction, has been recently 
discussed by comparing QMC calculations of the binding energy of hypernuclei 
with available experimental measurements (101 102). The role of A hyperons in 
neutron star matter has been discussed by Lonardoni et al. in Ref. (103). The 
main conclusion is that a small change of the A-neutron-neutron interaction, while 
producing very small changes to the properties of hypernuclei, has a dramatic 
effect in neutron matter. 

In Fig. m we show the mass-radius of neutron stars obtained by including 
the A hyperons to the EOS. As described in Ref. (103), the fraction of A is 
determined by imposing chemical equilibrium once the energy as a function of 
the concentrations of A is known. That has been calculated using the AFDMC 
method for different Hamiltonians. The AV^+UIX for the pure neutron sector 
has been combined with the A-nucleon interaction of Ref. (104), and with the 
addition of two different A-nucleon-nucleon potentials, denoted as AN+ANN (I) 
and AN+ANN (II). While the operatorial structure of these forces is the same, 
the strengths have been slightly changed to reproduce better the binding energy 
of A-hypernuclei. Although the two Hamiltonians give qualitatively very similar 
results in hypernuclei, the EOS and the corresponding neutron star structure are 
dramatically different. It is clear that further input on hyperon-nucleon interac¬ 
tions, especially on ANN, will be crucial if firm conclusions are to be drawn. 


6 CONCLUSIONS AND OUTLOOK 

A great deal has been learned about neutron matter over the past decade, both 
from theory and through experiments and observations. The equation of state 
of cold neutron matter at low and moderate densities is now tightly constrained 
from many-body calculations with realistic interactions. Also, the pairing gap 
has been reliably extracted from methods that correctly predicted the pairing 
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Figure 11: The mass-radius for neutron stars obtained from EOS that include A- 
hyperons. The green upper line (denoted with PNM) is the result of pure neutron 
matter (the same presented in Fig. 10). The red bottom band (indicated with 
AN) is obtained by including only the A-nucleon interaction in the Hamiltonian. 
The blue curve, AN+ANN (I), is obtained but including the Usmani three-body 
interaction (104). The EOS obtained using the AN+ANN (II) does not allow to 
A hyperons to produce, and is indicated by the dashed black line. Recent heavy- 
mass neutron stars measurements are also indicated by horizontal lines. Figure 
taken from (103). 


gap in Fermi atoms at unitarity, a similarly strongly paired Fermi system. 

Calculations of the equation of state provide strong constraints on the radii 
of typical neutron stars of masses ~ 1.4 solar masses. These seem to be in 
reasonable agreement with mass/radius constraints extracted from astrophysical 
observations, but the observations and particularly the associated extraction of 
the mass/radius relation remain a controversial area with a variety of results. 

Constraints on inhomogeneous matter are also starting to appear: it is encour¬ 
aging that these results can be made consistent with density functionals that also 
describe nuclei. Traditionally, density functionals have used limited constraints 
at the extremes of isospin, so this remains a promising avenue for future research. 
It will be interesting to see if these constraints alter predictions for the properties 
of neutron-star crusts. 

Finally, the observation of two-solar-mass neutron stars has challenged the 
traditional view that hyperons and other non-nucleonic degrees of freedom will 
appear and considerably soften the high-density equation of state. In the case 
of hyperons, the role of the hyperon-two-nucleon interaction is very important, 
and may limit the role of hyperons in dense neutron-rich matter. Studies of 
cold dense matter from the fundamental degrees of freedom, quarks and gluons, 
remain critical, though they are quite challenging. 
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There are many prospects for important further advances, both in theory and 
in experiments and observations. Experiments at rare isotope facilities will allow 
us to study many important neutron-rich nuclei, providing valuable information 
on neutron-rich matter. Observations and improved understanding should enable 
better constraints on the neutron-star mass-radius relation. The observation of 
gravitational waves associated with neutron-star mergers is an exciting possibility, 
and this could provide important constraints on neutron-star properties in the 
near future. 

Our increasing ability to calculate properties of strongly interacting quantum 
systems will also play an important role. Neutron-star matter typically contains 
approximately ten percent protons: the impact on the neutron star mass/radius 
relation is likely to be modest but it may be very important for dynamic properties 
including the weak response of dense neutron matter. The role of superfluidity at 
higher densities, both p-wave superfluidity in neutrons and s-wave superfluidity 
of the low-density protons, is another important problem. Finally, studies of 
matter at finite temperature are critical for both potential sites of heavy-element 
synthesis, core-collapse supernovae, and neutron-star mergers. 
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